/*
	Velocity2: A class that generates a rigib body rotation within the grid
			   It is used for testing the functionality of the library
*/

#ifndef VELOCITY_H
#define VELOCITY_H

//#include "main.h"
#include "vector.h"
#include "grid.h"
#include "matrix.h"
class Velocity
{
public:
   int cont;
   Grid gridU;
   Grid gridV;
   Grid gridW;
  // Matrix <double> U;
   //Matrix <double> V;
   //Matrix <double> W;
	Velocity(int nx,int ny, int nz):gridU(nx, ny, nz),gridV(nx, ny, nz),gridW(nx, ny, nz)
	{	}
   
	inline void GetVelocity(const Vector &pos, Vector &u) const
	{
		//vortex around 0,0,1
		u = Vector( gridU[int (pos[0])],gridV[int (pos[1])], gridW[int (pos[2])]);
      //if(cont ==0) cout<< "(" << u[0] << "," << u[1] << "," << u[2] << ")  ";
	}
};
class Velocity_Vector
{
    public:
    typedef double* pdouble;
    pdouble gridUVW[3];
    int Nx4,Nx4Ny4;
    int size_comp,size_comp_x_n_copies;
    Velocity_Vector(int nx,int ny, int nz, int n_copies):Nx4(nx+4),Nx4Ny4((nx+4)*(ny+4))
    { 
	 size_comp = (nx+4) * (ny+4) * (nz+4);
         size_comp_x_n_copies  = n_copies *  size_comp;
         
	 for(int i=0; i<=2; i++)
         {
             gridUVW[i] = new double[size_comp_x_n_copies]; 
             fill(gridUVW[i], gridUVW[i]+size_comp_x_n_copies,0.0);
	     
         }
         
    }
        inline Vector operator[](int index){return Vector(*(gridUVW[0]+index),*(gridUVW[1]+index),*(gridUVW[2]+index));}
	inline double& operator()(int h, int index){return *(gridUVW[h]+index);}
	inline Vector operator()(int i, int j, int k){int index = i+j*Nx4+k*(Nx4Ny4); return (*this)[index];}
	inline double& operator()(int h, int i, int j, int k) {int index = i+j*Nx4+k*(Nx4Ny4);return (*this)(h,index);}
	inline void GetVelocity_Vector(const Vector &pos, Vector &u) 
	{
		u=(*this)(int (pos[0]),int (pos[1]),int (pos[2]));
	}
};

#endif
